Fonte: Storymap UFRRJ

Introdução à Análise Estatística Espacial

O que é Análise Estatística Espacial ?

  • São métodos estatísticos que levam em consideração a localização espacial do fenômeno estudado;

  • Segundo Bailey & Gatrell (1995), “Análise estatística espacial pode ser realizada quando os dados são espacialmente localizados e se considera explicitamente a possível importância de seu arranjo espacial na análise ou interpretação dos resultados”

  • A primeira pergunta a ser feita é: A distribuição dos dados apresenta um padrão aleatório ou apresenta uma agregação definida (clusters) ?

Origem da Estatística Espacial

  • Dr. John Snow (1813-1858) \(\rightarrow\) Considerado pai da Epidemiologia Moderna

Mapeamento dos casos de coléra (\(\bullet\)) e as bombas de água (X) em Londres, 1854.

Objetivos da Estatística Espacial

  1. Análise de padrões espaciais e espaço-temporais (ex: Análise Exploratória Espacial, Correlação Espacial)

  2. Modelagem Espacial (ex: Modelos Estatísticos de Regressão Espacial)

Dependência Espacial ou Autocorrelação Espacial

  • “Independência é um pressuposto muito conveniente que faz grande parte da teoria estatística matemática tratável. Entretanto, modelos que envolvem dependência estatística são frequentemente mais realísticos. . . . Neste tipos de dados a dependência está presente em todas as direções e fica mais fraca a medida em que aumenta a dispersão na localização dos dados.” (Cressie N.,1991)

  • “Todas as coisas se parecem, coisas mais próximas são mais parecidas que aquelas mais distantes.” (Tobler, 1979). Também conhecida como \(1^a\) Lei da Geografia

Tipologia dos Dados Espaciais

Os diferentes tipos de dados espaciais são tradicionalmente classificados de acordo com uma tipologia. Esta caracterização diz respeito a natureza estocástica da observação.

Noel Cressie (1993) divide a estatı́stica espacial em 3 grandes áreas:

  1. Dados de processos pontuais

  2. Dados de geoestatística

  3. Dados de área

Existem métodos estatı́sticos diferentes para descrever ou analisar estes tipos de dados. Exemplo:

Padrões Pontuais

  • O principal interesse está no conjunto de coordenadas geográficas representando as localizações exatas de eventos.

  • Exemplos: Localização de crimes, localização da residência dos casos de dengue, localização de espécies vegetais, etc.

  • Neste caso, o dado aleatório de interesse é a localização espacial do evento.

Estimativas de Kernel (ou Mapas de Calor)

  • Estimador de intensidade de distribuição de pontos:

\[\hat{\lambda}_{\tau}(u) = \dfrac{1}{\tau^2}\sum k(\dfrac{d(u_i , u)}{\tau}), d(u_i , u) \leq \tau\] *Fonte: Referência científica: Druck, S.; Carvalho, M.S.; Câmara, G.; Monteiro, A.V.M. (eds) “Análise Espacial de Dados Geográficos”. Brasília, EMBRAPA, 2004

  • Localização da ocorrência de casos de dengue em Belo Horizonte/MG

Geoestatística

  • Podemos definir como sendo uma análise de um atributo espacialmente contínuo amostrado em localizações fixas.

  • Os dados compreendem um conjunto de localizações (em geral latitudes e longitudes), mas atribuidos a eles uma medida, como por exemplo: volume de chuva, concentração de poluentes no ar, número de ovos de Aedes postado em ovitrampas, etc.

  • A determinação experimental do semivariograma, para cada valor de \(x_i\), considera todos os pares de amostras \(z(x)\) e \(z(x+h)\), separadas pelo vetor distância \(h\), a partir da equação:

\[\hat{\gamma}(h) = \dfrac{1}{2N(h)}\sum^{N(h)}_{i=1}[z(x_i) - z(x_i + h)]^2 \]

  • Sendo \(\hat{\gamma}(h)\) o semivariograma estimado e \(N(h)\) o número de pares de valores medidos, \(z(x)\) e \(z(x+h)\), separados pelo vetor \(h\).

  • Esta fórmula, entretanto, não é robusta. Podem existir situações em que variabilidade local não é constante e se modifica ao longo da área de estudo (heteroscedasticidade).

*Fonte: Referência científica: Druck, S.; Carvalho, M.S.; Câmara, G.; Monteiro, A.V.M. (eds) “Análise Espacial de Dados Geográficos”. Brasília, EMBRAPA, 2004

  • Exemplo: Mapa sobre o teor de argila no solo.

Dados de Área

  • Na análise de áreas o atributo estudado é em geral resultado de uma medida sumáraria (ex: contagem, taxas, médias, etc.) em uma determinada área bem definida, por exemplo um polígono.

  • Tal polígono cuja forma pode ser complexa bem como as relações de vizinhança.

  • O objetivo é a detecção e explicação de padrões e tendências observados entre áreas.

Mapa Temático

  • Taxas de câncer de pulmão na população branca masculina nos Estados Unidos, por condados no ano de 1998

Matriz de Vizinhança

Moran Global, Moran Local e Lisa Map

Moran Global

Moran Global

Moran Local

Moran Local

  • Desigualdade no nível distrital na cobertura de saúde reprodutiva, materna, neonatal e infantil na Índia

Fonte: PANDA, Basant Kumar; KUMAR, Gulshan; AWASTHI, Ashish. District level inequality in reproductive, maternal, neonatal and child health coverage in India. BMC public health, v. 20, n. 1, p. 1-10, 2020.

Modelos de Regressão Espacial

  • Hipótese de independência das observações em geral é Falsa \(\rightarrow\) Dependência Espacial

  • Efeitos Espaciais \(\rightarrow\) Se existir forte tendência ou correlação espacial, os resultados serão influenciados, apresentando associação estatística onde não existe (e vice-versa);

  • Como verificar ? \(\rightarrow\) Medir a autocorrelação espacial dos resíduos da regressão (Índice de Moran dos resíduos).

  • Autocorrelação espacial constatada ! E agora ? \(\rightarrow\) Utilizar Modelos de regressão que incorporam efeitos espaciais:

    1. Modelos Globais: utilizam um único parâmetro para capturar a estrutura de correlação espacial \(\rightarrow\) ex: Spatial Error Models (CAR)

\[y_i = \beta_0 + \sum^{p}_{k} \beta_k x_{ik} + \varepsilon_i \text{ , sendo } \varepsilon_i = \lambda W + \xi\]

  • Os efeitos da autocorrelação espacial são associados ao termo de erro. Sendo \(W\) a matriz de vizinhaça, \(\varepsilon\) a componente do erro com efeitos espaciais, \(λ\) é o coeficiente autoregressivo e \(\xi\) é a componente do erro com variância constante e não correlacionada.

  • A hipótese nula para a não existência de autocorrelação é que \(\lambda = 0\), ou seja, o termo de erro não é espacialmente correlacionado.

    1. Modelos Locais: parâmetros variam continuamente no espaço ex: Geographically Weighted Regression (GWR)

\[y_i = \beta_{0}(u_i, v_i) + \sum^{p}_{k} \beta_k (u_i, v_i) x_{ik} + \varepsilon_i \text{ , sendo } (u_i, v_i) \text{ as coordenadas geográficas}\] Fonte: ARDILLY, Pascal et al. Manuel d’analyse spatiale. 2018.

Serão feitas quatro aplicações da Estatística Espacial utilização o pacote estatístico R:

Aplicação I: Utilizando a biblioteca tmap para construção de mapas temáticos

library(tmap)
data(World)
tmap_style("classic")
# Desenhando um rápido mapa temático mundial para esperança de vida.
qtm(World, fill = "life_exp")

Aplicação II: Baixando e construindo mapas a partir da biblioteca geobr

  • Bibliotecas que serão utilizadas:
library(ggplot2)
library(dplyr)
library(viridis)
library(geobr)
library(sf)
library(maptools)
library(leaflet)
library(ggspatial)
  • Para acessar os dados dos limites territoriais de todos os estados brasileiros é necessário utilizar a função read_state.
brasil.ufs <- read_state(code_state = "all", year=2019, showProgress = FALSE)
  • Primeiro, vamos fazer um gráfico apenas com as geometrias.
ggplot(brasil.ufs) + 
  geom_sf()

  • Para a construção de um mapa onde cada estado recebe uma cor de acordo com a sua região geogrática, procedemos da seguinte forma:
ggplot(brasil.ufs) + 
  geom_sf(aes(fill = name_region)) + 
  labs(fill="Região")

  • Agora utilizaremos dados relativos a porcentagem de municípios com rede de esgoto de acordo com a unidade da federação (fonte dos dados ) para fazer um gráfico. Vamos associar esses dados a tabela de acordo com a variável State e padronizaremos a porcentagem para variar de 0 a 100.
# Entrando com os dados observados na wikipedia
acesso_san <- data.frame(code_state = c(12, 27, 16, 13, 29, 23, 53, 32, 52, 21, 51, 50, 31, 15, 
                                   25, 41, 26, 22, 33, 24, 43, 11, 14, 42, 35, 28, 17), 
                         com_rede = c(0.273, 0.412, 0.313, 0.177, 0.513, 0.696, 1.000, 0.974, 0.280, 0.065, 
                                      0.191, 0.449, 0.916, 0.063, 0.731, 0.421, 0.881, 0.045, 0.924, 0.353, 
                                      0.405, 0.096, 0.400, 0.352, 0.998, 0.347, 0.129))
# Construindo o mapa com ggplot
brasil.ufs |> 
  left_join(acesso_san, by = "code_state") |> 
  mutate(com_rede = 100*com_rede) |> 
  ggplot(aes(fill = com_rede), color = "black") +
    geom_sf() + 
    scale_fill_viridis(name = "Municípios com rede de esgoto (%)", direction = -1) + 
  xlab("Longitude") + 
  ylab("Latitude") +
  annotation_scale(location = "bl") +
  annotation_north_arrow(location = "br") +
    theme_minimal()

  • Uma forma alteranativa de apresentar esses mesmos dados se dá pela apresentação de círculos com raios proporcionais a porcentgem de municípios com rede de esgoto no centroide de cada geometria (nesse caso, UF).
# Juntando o banco com os dados + malha
coord_pontos <- brasil.ufs |> 
                  left_join(acesso_san, by = "code_state") |> 
                  mutate(com_rede = 100*com_rede) |> 
                  st_centroid()
# Construindo o mapa com ggplot
ggplot(brasil.ufs)+ 
  geom_sf() + 
  geom_sf(data = coord_pontos, aes(size = com_rede), col = "blue", alpha = .65,
          show.legend = "point") + 
  scale_size_continuous(name = "Municípios com rede de esgoto (%)") + 
  xlab("Longitude") + 
  ylab("Latitude") +
  annotation_scale(location = "bl") +
  annotation_north_arrow(location = "br") +
    theme_minimal()

  • Uma alternativa interativa para trabalhar com mapas é com a utilização do pacote leaflet
data.frame(st_coordinates(coord_pontos), 
           com_rede = coord_pontos$com_rede, 
           UF = coord_pontos$name_state) |>
  leaflet() |> 
    addTiles() |>
    addCircleMarkers(~ X, ~ Y,
                     label = ~ as.character(paste0(UF, ": ", com_rede, "%")),
                     labelOptions = labelOptions(textsize = "13px"),
                     radius = ~ com_rede/10,
                     fillOpacity = 0.5)

Créditos ao Thiago Mendonça

Aplicação III: Dengue em Dourados/MS - Parte 1: Análise exploratória

OBS: Os dados e as malhas geográficas utilizadas nessa apresentação, estão disponíveis no seguinte endereço: (link)

  • Biliotecas do R que serão utilizadas
library(readr)
library(tidyverse)
library(sf)
library(maptools)
library(spatstat)
library(tmap)
  • Lendo a tabela da população por setor censitário e os shapes files do contorno e por setor censitário de Dourados/MS
unzip('dados/pop2010.zip', exdir='dados')
pop2010 <- read_csv('dados/pop2010.csv')

unzip('malhas/Setor_UTM_SIRGAS.zip', exdir='malhas')
setor.sf <- read_sf('malhas/Setor_UTM_SIRGAS.shp', crs = 31981)

unzip('malhas/contorno.zip', exdir='malhas')
contorno.sf <- read_sf('malhas/contorno.shp', crs = 31981)

OBS: Um aspecto muito importante dos dados espaciais é o sistema de referência de coordenadas (CRS) que é usado. Por exemplo, uma localização de (140, 12) não é significativa se você sabe onde está a origem e se a coordenada x está a 140 metros, quilômetros ou talvez graus de distância dela (na direção x). (link)

  • Fazendo um join com as tabelas com os setores censitários + população
setor.sf <- setor.sf %>% mutate (idsetor = as.numeric(CD_GEOCODI)) %>% inner_join(pop2010,by='idsetor')
  • Lendo e plotando os casos de dengue georreferenciados em Dourados/MS
unzip('dados/dengue_dourados.zip', exdir='dados')
casos <- read_csv('dados/dengue_dourados.csv')
casos.sf <- st_as_sf(casos, coords = c("X", "Y"), crs = 31981)
ggplot(setor.sf) + 
  geom_sf(fill = 'white', color='black') +
  geom_sf(data=casos.sf, color='red',size=1) +
  theme_void()

  • Lendo e plotando os os pontos de coleta de lixo georreferenciados em Dourados/MS
unzip('dados/lixo_dourados.zip', exdir='dados')
lixo <- read_csv2('dados/lixo_dourados.csv')
lixo.sf <- st_as_sf(lixo, coords = c("X", "Y"), crs = 31981)
ggplot(setor.sf) + 
  geom_sf(fill = 'white', color = 'black') +
  geom_sf(data=lixo.sf,color='blue',size=1) +
  theme_void()

Como podemos observar, existem alguns pontos de coleta de lixo fora do contorno de Dourados/MS

  • Uma forma de ficarmos só com os pontos dentro do polígono é utilizando o comando st_intersection.
lixo2.sf <- st_intersection(contorno.sf, lixo.sf)
ggplot(setor.sf) + 
  geom_sf(fill = 'white', color = 'black') +
  geom_sf(data=lixo2.sf,color='blue',size=1) +
  theme_void()

  • Utilizando as informações dos casos (pontos) + do lixo (ponto) + população de cada setor censitário (mapa temático)
ggplot(setor.sf) + 
  geom_sf(aes(fill=pop)) + 
  geom_sf(data=casos.sf,color='red',size=0.7, aes(colour = "Caso"), 
          show.legend = "point") +
  geom_sf(data=lixo2.sf,color='salmon',size=1, aes(colour = "Lixo"), 
          show.legend = "point") +
  scale_fill_distiller(palette ="PuBu", direction = 1) +
  scale_colour_manual(values = c("Caso" = "red", "Lixo" = "slamon")) +
  theme_minimal()

  • Agora serão construidos buffers de 500m de distância ao redor de cada ponto de coleta de lixo. Isso é interessante para verificar se os casos de dengue ocorrem em um raio de até 500m de distância dos pontos de coleta de lixo. Ou seja, a pergunta é, será que a distância dos pontos de coeta de lixo influenciam a ocorrência do caso de dengue ?

Buffers: São polígonos que contornam um objeto a uma determinada distância. Sua principal função é materializar os conceitos de “perto” e “longe”.

lixo_buffer <- st_buffer(lixo2.sf, 500) 

ggplot(setor.sf) + 
  geom_sf(aes(fill=pop)) + 
  geom_sf(data=lixo_buffer,color='gray', fill = "transparent", size=0.4) +
  geom_sf(data=casos.sf, color='red',size=0.7) +
  scale_fill_distiller(palette ="PuBu", direction = 1) +
  scale_colour_manual(values = c("Caso" = "red", "Lixo" = "slamon")) +
  theme_minimal()

  • Represntando os casos e o lixo de forma interativa
tm_shape(setor.sf) +  tm_borders("black") +
  tm_shape(casos.sf) + tm_dots("red") +
  tm_shape(lixo.sf) + tm_dots("green") +
  tm_shape(lixo_buffer) + tm_borders("blue") +
  tmap_mode("view")
  • Convertendo o dado de pontos (padrão pontual) para dados de área
## conta casos por setor 
casos.sf$contador <- 1 
casos <- setor.sf |>  
  st_join(casos.sf) |> 
  filter(CLASSI_FIN == 1) %>%  ## seleciona somente os casos confirmados
  group_by(ID) |> 
  summarise(casos = sum(contador))

st_geometry(casos) <- NULL  ## remove atributos de geometria

## numero de depositos de Lixo por setor 
lixo.sf$contador <- 1 

lixo <- setor.sf |>  
  st_join(lixo.sf) |> 
  group_by(ID) |> 
  summarise(lixo = sum(contador))

st_geometry(lixo) <- NULL ## remove atributos de geometria

# Inserindo as contagens dos casos e de pontos de coleta de lixo no atributo com a geometria.
setor.sf <- setor.sf |> 
  left_join(lixo,by = 'ID') |> 
  left_join(casos,by = 'ID')
  • Plotando o mapa temático dos casos por setor censitário
plot(setor.sf['casos'])

  • Plotando o mapa temático dos pontos de coleta de lixo por setor censitário
plot(setor.sf['lixo'])

  • Calculando a taxa de incidência e plotando o mapa temático dos pontos de coleta de lixo por setor censitário
setor.sf$tx <- (setor.sf$casos/setor.sf$pop) * 1000
setor.sf$tx[is.na(setor.sf$tx)] <- 0 # Transformando os missings em zero

summary(setor.sf$tx)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   0.000   1.736   4.065   4.311  56.399
  • Plotando a distribuição da incidência em Dourados/MS
library(wesanderson)
pal <- wes_palette("Moonrise3", 20, type = "continuous")

ggplot(setor.sf) + 
  geom_sf(aes(fill = tx), color = 'black') +
  scale_fill_gradientn(colours = pal) +
  ggtitle("Taxa de incidência de Dengue") + 
  theme_void()

  • Kernel por atributos

Vamos plotar o kernel por atributos referente a taxa de incidência de dengue em Dourados/MS.

Primeiramente é necessário dissolver os polígonos em formato sf para obter o contorno. Nesse caso queremos preservar o atributo AREA

dourados.contorno <- st_union(setor.sf)
plot(dourados.contorno)

dourados.w <- as.owin(st_geometry(dourados.contorno))
  • Extraindo os centróides dos polígonos em Dourados/MS
centroides <- st_centroid(st_geometry(setor.sf))

# Transformando em os centróides em formato sp
centroides.sp <- as.data.frame(as_Spatial(centroides))
names(centroides.sp) <- c('X','Y')

plot(centroides.sp)

  • Colocando os pontos no formato sp
centroides.ppp <- ppp(centroides.sp$X,centroides.sp$Y, dourados.w)

plot(centroides.ppp,pch=19,cex=0.5)

  • Fazendo o kernel por atributo da taxa de detecção
kernel.tx <- density(centroides.ppp, 500, weights = setor.sf$tx, scalekernel = TRUE)
plot(kernel.tx)

  • Construindo a matriz de vizinhança para verificar a autocorrelação espacial
library(spdep)
viz <- poly2nb(setor.sf)
viz 
Neighbour list object:
Number of regions: 284 
Number of nonzero links: 1726 
Percentage nonzero weights: 2.139952 
Average number of links: 6.077465
  • Iremos precisar da coordenadas dos centróides
setor.sp <- as(setor.sf, 'Spatial') # convertendo em formato sp
coord <- coordinates(setor.sp) # coordenadas dos centroidas dos poligonos de dourados
class(setor.sp)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
  • Verificando a malha de conectividade da vizinhança de Dourados/MS
viz.sf <- as(nb2lines(viz, coords = coord), 'sf')
viz.sf <- st_set_crs(viz.sf, st_crs(setor.sf))

# Plota o grafo de conectividade por contiguidade
mapa.viz <- ggplot(setor.sf) + 
  geom_sf(fill = 'salmon', color = 'white') +
  geom_sf(data = viz.sf) +
  theme_minimal() +
  ggtitle("Vizinhança por \n conectividade") +
  ylab("Latitude") +
  xlab("Longitude")
mapa.viz

  • Obtendo a correlação da taxa de incidência de dengue Dourados/MS
pesos.viz <- nb2listw(viz)
moran.test(setor.sf$tx, pesos.viz)

    Moran I test under randomisation

data:  setor.sf$tx  
weights: pesos.viz    

Moran I statistic standard deviate = 15.724, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.524720303      -0.003533569       0.001128660
  • Plotando o correlograma
correl <- sp.correlogram(viz, setor.sf$tx, order = 8, method = "I")
correl
Spatial correlogram for setor.sf$tx 
method: Moran's I
           estimate expectation    variance standard deviate Pr(I) two sided
1 (284)  0.52472030 -0.00353357  0.00112866          15.7239       < 2.2e-16
2 (284)  0.23776816 -0.00353357  0.00048540          10.9525       < 2.2e-16
3 (284)  0.09536593 -0.00353357  0.00028358           5.8729       4.282e-09
4 (284) -0.02063378 -0.00353357  0.00019736          -1.2172         0.22351
5 (284) -0.07589158 -0.00353357  0.00016732          -5.5939       2.221e-08
6 (284) -0.06448448 -0.00353357  0.00016165          -4.7940       1.635e-06
7 (284) -0.02708340 -0.00353357  0.00016725          -1.8210         0.06861
8 (284) -0.02744543 -0.00353357  0.00018328          -1.7662         0.07735
           
1 (284) ***
2 (284) ***
3 (284) ***
4 (284)    
5 (284) ***
6 (284) ***
7 (284) .  
8 (284) .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(correl)

  • Mapeando os polígonos que tiveram os p-valores mais significativos no Moran Local.
setor.sf$pval <- localmoran(setor.sf$tx, pesos.viz)[,5]

tm_shape(setor.sf) + 
  tm_polygons(col='pval', title="p-valores", breaks=c(0, 0.01, 0.05, 0.10, 1), border.col = "white", palette="-Oranges") +
  tm_scale_bar(width = 0.15)
  • Moran Local (Lisa Map) da taxa de incidência Dourados/MS
resI <- localmoran.sad(lm(setor.sf$tx ~ 1), 1:length(viz), viz, style = "W")
summary(resI)[1:10,]
    Local Morans I Stand. dev. (N)   Pr. (N) Saddlepoint Pr. (Sad)  Expectation
1 1   -0.009885292     -0.01575255 1.0125682 -0.03600714 0.9712767 -0.003533569
2 2    0.226981101      0.46511580 0.6418485  0.70056077 0.4835772 -0.003533569
3 3   -0.010910944     -0.01829621 1.0145974 -0.04041673 0.9677609 -0.003533569
4 4    0.484675519      1.31014141 0.1901480  1.43930439 0.1500643 -0.003533569
5 5    0.018648578      0.05952726 0.9525322  0.09384037 0.9252360 -0.003533569
     Variance    Skewness Kurtosis   Minimum  Maximum         omega       sad.r
1 1 0.1625854 -0.05147857 8.753481 -57.75455 56.75455 -0.0001369880 -0.01569409
2 2 0.2456263 -0.04188294 8.752890 -70.87400 69.87400  0.0028119325  0.44472643
3 3 0.1625854 -0.05147857 8.753481 -57.75455 56.75455 -0.0001590844 -0.01822753
4 4 0.1388594 -0.05570264 8.753780 -53.41199 52.41199  0.0066304485  1.08297547
5 5 0.1388594 -0.05570264 8.753780 -53.41199 52.41199  0.0005595065  0.05929966
          sad.u
1 1 -0.01569909
2 2  0.49831660
3 3 -0.01823490
4 4  1.59298211
5 5  0.05942125
 [ reached 'max' / getOption("max.print") -- omitted 5 rows ]
setor.sf$MoranLocal <- summary(resI)[,1] 

library(scales)

ggplot(setor.sf) + 
  geom_sf(aes(fill = MoranLocal), color = 'black') +
  scale_fill_gradientn(colours=c("blue", "white", "red"), 
                       values=rescale(c(min(setor.sf$MoranLocal), 0, max(setor.sf$MoranLocal))), guide="colorbar") + 
  ggtitle("Moran local") + 
  theme_void()

Aplicação IV: Dengue em Dourados/MS - Parte 2: Modelagem (Modelos Linear, CAR e GWR)

  • Ajustando o modelo de regressão linear simples.
# Transformando os missings em zero
setor.sf$lixo[is.na(setor.sf$lixo)] <- 0  

dourados.lm <- lm(tx ~ lixo, data = setor.sf)
summary(dourados.lm)

Call:
lm(formula = tx ~ lixo, data = setor.sf)

Residuals:
   Min     1Q Median     3Q    Max 
-4.510 -4.115 -2.286  0.237 51.889 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   4.5103     0.4889   9.225   <2e-16 ***
lixo         -0.3955     0.1945  -2.033    0.043 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.365 on 282 degrees of freedom
Multiple R-squared:  0.01445,   Adjusted R-squared:  0.01095 
F-statistic: 4.134 on 1 and 282 DF,  p-value: 0.04298
  • Checando os residuos para verificar a presença de autocorrelação.
dourados.lm$lmresid <- residuals(dourados.lm)
moran.test(dourados.lm$lmresid, pesos.viz)

    Moran I test under randomisation

data:  dourados.lm$lmresid  
weights: pesos.viz    

Moran I statistic standard deviate = 15.578, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.520064269      -0.003533569       0.001129669
  • Ajustando o modelo CAR (Spatial Error Model)
library(spatialreg)
dourados.car <- errorsarlm(tx ~ lixo, data = setor.sf, listw = pesos.viz)
summary(dourados.car)

Call:errorsarlm(formula = tx ~ lixo, data = setor.sf, listw = pesos.viz)

Residuals:
      Min        1Q    Median        3Q       Max 
-12.53213  -2.34032  -1.09686   0.82021  31.62992 

Type: error 
Coefficients: (asymptotic standard errors) 
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  4.07199    1.54998  2.6271 0.008611
lixo        -0.25129    0.14814 -1.6963 0.089826

Lambda: 0.8063, LR test value: 167.98, p-value: < 2.22e-16
Asymptotic standard error: 0.041575
    z-value: 19.394, p-value: < 2.22e-16
Wald statistic: 376.13, p-value: < 2.22e-16

Log likelihood: -885.0643 for error model
ML residual variance (sigma squared): 25.435, (sigma: 5.0433)
Number of observations: 284 
Number of parameters estimated: 4 
AIC: NA (not available for weighted model), (AIC for lm: 1944.1)
  • Checando os residuos para verificar a presença de autocorrelação
dourados.car$carresid <- residuals(dourados.car)
moran.test(dourados.car$carresid, pesos.viz)

    Moran I test under randomisation

data:  dourados.car$carresid  
weights: pesos.viz    

Moran I statistic standard deviate = 0.78357, p-value = 0.2166
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.023005079      -0.003533569       0.001147099
  • Ajustando o modelo GWR (Geographically Weighted Regression)
# Biblioteca para ajustar o modelos GWR
library(spgwr)
# Estimando a largura de banda “ideal” para o kernel
GWRbanda <- gwr.sel(tx ~ lixo, data = setor.sf, coords = cbind(centroides.sp$X, centroides.sp$Y), 
    adapt = T)
Adaptive q: 0.381966 CV score: 14699.46 
Adaptive q: 0.618034 CV score: 15071.4 
Adaptive q: 0.236068 CV score: 14204.85 
Adaptive q: 0.145898 CV score: 13336.85 
Adaptive q: 0.09016994 CV score: 12151.44 
Adaptive q: 0.05572809 CV score: 10995.75 
Adaptive q: 0.03444185 CV score: 10093.37 
Adaptive q: 0.02128624 CV score: 9758.837 
Adaptive q: 0.01315562 CV score: 9295.706 
Adaptive q: 0.008130619 CV score: 8996.895 
Adaptive q: 0.005024999 CV score: 8988.675 
Adaptive q: 0.00638842 CV score: 9000.179 
Adaptive q: 0.00310562 CV score: 9842.835 
Adaptive q: 0.004291861 CV score: 9067.819 
Adaptive q: 0.005545779 CV score: 8976.568 
Adaptive q: 0.005867639 CV score: 8980.638 
Adaptive q: 0.005586469 CV score: 8976.67 
Adaptive q: 0.005505089 CV score: 8976.598 
Adaptive q: 0.005545779 CV score: 8976.568 
# Ajustando o modelo GWR
dourados.gwr = gwr(tx ~ lixo, data = setor.sf, coords = cbind(centroides.sp$X, centroides.sp$Y), 
    adapt = GWRbanda, hatmatrix = TRUE, se.fit = TRUE)

dourados.gwr
Call:
gwr(formula = tx ~ lixo, data = setor.sf, coords = cbind(centroides.sp$X, 
    centroides.sp$Y), adapt = GWRbanda, hatmatrix = TRUE, se.fit = TRUE)
Kernel function: gwr.Gauss 
Adaptive quantile: 0.005545779 (about 1 of 284 data points)
Summary of GWR coefficient estimates at data points:
                   Min.    1st Qu.     Median    3rd Qu.       Max.  Global
X.Intercept.  -0.324952   1.412188   3.099163   5.397652  37.581278  4.5103
lixo         -18.321709  -1.222217  -0.400929  -0.061492   1.555744 -0.3955
Number of data points: 284 
Effective number of parameters (residual: 2traceS - traceS'S): 134.4484 
Effective degrees of freedom (residual: 2traceS - traceS'S): 149.5516 
Sigma (residual: 2traceS - traceS'S): 5.447405 
Effective number of parameters (model: traceS): 100.4864 
Effective degrees of freedom (model: traceS): 183.5136 
Sigma (model: traceS): 4.917575 
Sigma (ML): 3.952993 
AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 1904.233 
AIC (GWR p. 96, eq. 4.22): 1687.144 
Residual sum of squares: 4437.827 
Quasi-global R2: 0.7140881 
# Colocando a saída do modelo dentro de um objeto dataframe.
results <- as.data.frame(dourados.gwr$SDF)
head(results)
     sum.w X.Intercept.       lixo X.Intercept._se  lixo_se       gwr.e
1 4.630090     6.489949 -1.8954895        2.411594 1.718037 -0.02159391
2 3.934954     7.883872 -1.2811203        2.504116 1.630176  1.98454873
3 3.062602     7.267798 -1.5625806        2.780695 1.625888 -1.85906371
4 4.139484     9.791258 -1.1687314        2.376803 1.277302  6.80957181
5 5.541437     4.586352 -0.6839813        2.090863 1.185501 -2.63010618
      pred  pred.se   localR2 X.Intercept._se_EDF lixo_se_EDF pred.se.1
1 2.698970 2.477284 0.6403656            2.671425    1.903141  2.744191
2 7.883872 2.504116 0.5043081            2.773914    1.805815  2.773914
3 5.705218 2.164469 0.5288378            3.080292    1.801064  2.397673
4 8.622527 1.794393 0.4101644            2.632885    1.414921  1.987725
5 3.902371 1.524868 0.6665163            2.316137    1.313229  1.689160
   coord.x coord.y
1 731406.5 7541547
2 730845.4 7541481
3 730957.7 7541239
4 730487.2 7541437
5 729874.1 7541446
 [ reached 'max' / getOption("max.print") -- omitted 1 rows ]
  • Verificando a distribuição dos coeficientes de regressão para a variável lixo
hist(results$lixo)
abline(v = median(results$lixo), col = "red")

  • Verificando a distribuição dos localR2
hist(results$localR2)
abline(v = median(results$localR2), col = "blue")

  • Incorporando alguns parâmetros de saída do modelo na tabela setor.sf
setor.sf$coef.lixo <- results$lixo
setor.sf$localR2 <- results$localR2
setor.sf$pred.gwr <- results$pred
  • Mapa dos coeficientes de regressão para a variável lixo
map.lixo <- ggplot(setor.sf) + 
            geom_sf(aes(fill = coef.lixo), color = "black") +              scale_fill_gradientn(colours = pal) + ggtitle("Distribuição dos coef. var. lixo") + 
            theme_void()
map.lixo

  • Checando os residuos para verificar a presença de autocorrelação para o modelo GWR.
# Calculando os resíduos para o modelo GWR
results$residuos <- setor.sf$tx - results$pred

moran.test(results$residuos, pesos.viz)

    Moran I test under randomisation

data:  results$residuos  
weights: pesos.viz    

Moran I statistic standard deviate = 1.4806, p-value = 0.06935
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.046816835      -0.003533569       0.001156432
  • Mapeando os coeficientes de regressão para a variável lixo por significancia através do teste de wald
# Calculando a estatística de wald
setor.sf$wald.teste <- abs(results$lixo/results$lixo_se)
# Dicotomizando a estatística de wald
setor.sf$wald.teste <- ifelse(setor.sf$wald.teste < 2, 0, 1)


map.wald <- ggplot(setor.sf) + geom_sf(aes(fill = factor(wald.teste)), color = "black") + 
    scale_fill_manual(values = c("white", "purple"), labels = c("< 2", ">= 2"), name = "Wald") + 
    ggtitle("Coef. lixo significativos") + theme_void()


library(gridExtra)
grid.arrange(map.lixo, map.wald, ncol = 2)

  • Mapa dos coeficientes de determinação regionalizados (R2 local).
ggplot(setor.sf) + 
  geom_sf(aes(fill = localR2), color = "black") + scale_fill_gradientn(colours = pal) + 
  ggtitle("R² local") + theme_void()

  • Verificando a distribuição dos preditos.
histdens <- function(x, titulo = "") {
    densi <- density(x)
    xli <- range(densi$x)
    yli <- range(densi$y)
    hist(x, col = "red", probability = T, xlim = xli, ylim = yli, main = titulo)
    lines(densi, lwd = 2)
    abline(v = median(x), lwd = 4, col = 4, lty = 2)
}

attach(mtcars)
par(mfrow = c(2, 2))

hist.tx <- histdens(setor.sf$tx, "Tx Bruta")
hist.lm <- histdens(dourados.lm$fitted.values, "Pred LM")
hist.car <- histdens(dourados.car$fitted.values, "Pred CAR")
hist.gwr <- histdens(results$pred, "Pred GWR")

  • Mapeando os valores observados e preditos dos modelos ajustados
library(colorspace)  # 

setor.sf$brks <- cut(setor.sf$tx, include.lowest = TRUE, right = TRUE, breaks = c(-4, 
    0, 2, 4, 10, 57), labels = c("0", "0 - 2", "2 - 4", "4 - 10", "> 10"))

tx.bruta <- ggplot(setor.sf) + geom_sf(aes(fill = brks), color = "black") + ggtitle("Taxa Bruta") + 
    scale_fill_discrete_sequential(palette = "Heat2", c1 = 80, c2 = 30, l1 = 30, 
        l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75", drop = FALSE, name = "Taxa") + 
    theme_void()



setor.sf$brks.lm <- cut(dourados.lm$fitted.values, lowest = TRUE, right = TRUE, breaks = c(-4, 
    0, 2, 4, 10, 57), labels = c("0", "0 - 2", "2 - 4", "4 - 10", "> 10"))


pred.lm <- ggplot(setor.sf) + geom_sf(aes(fill = brks.lm), color = "black") + ggtitle("Taxa Predita - LM") + 
    scale_fill_discrete_sequential(palette = "Heat2", c1 = 80, c2 = 30, l1 = 30, 
        l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75", drop = FALSE, name = "Taxa") + 
    theme_void()

setor.sf$brks.car <- cut(dourados.car$fitted.values, lowest = TRUE, right = TRUE, 
    breaks = c(-4, 0, 2, 4, 10, 57), labels = c("0", "0 - 2", "2 - 4", "4 - 10", 
        "> 10"))


pred.car <- ggplot(setor.sf) + geom_sf(aes(fill = brks.car), color = "black") + ggtitle("Taxa Predita - CAR") + 
    scale_fill_discrete_sequential(palette = "Heat2", c1 = 80, c2 = 30, l1 = 30, 
        l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75", drop = FALSE, name = "Taxa") + 
    theme_void()

setor.sf$brks.gwr <- cut(results$pred, lowest = TRUE, right = TRUE, breaks = c(-4, 
    0, 2, 4, 10, 57), labels = c("0", "0 - 2", "2 - 4", "4 - 10", "> 10"))


pred.gwr <- ggplot(setor.sf) + geom_sf(aes(fill = brks.car), color = "black") + ggtitle("Taxa Predita - GWR") + 
    scale_fill_discrete_sequential(palette = "Heat2", c1 = 80, c2 = 30, l1 = 30, 
        l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75", drop = FALSE, name = "Taxa") + 
    theme_void()

library(gridExtra)
grid.arrange(tx.bruta, pred.lm, pred.car, pred.gwr, ncol = 2)

  • Verificando a distribuição dos resíduos
library(vioplot)
vioplot(dourados.lm$residuals, dourados.car$residuals, results$residuos, names = c("LM", 
    "CAR", "GWR"), col = "orange")
title("Gráficos violinos da distribuição dos resíduos")

  • Mapeando os resíduos dos modelos ajustados
library(colorspace)  # 

setor.sf$brks.res.lm <- cut(dourados.lm$residuals, include.lowest = TRUE, right = TRUE, 
    breaks = c(-14, -5, -1, 1, 5, 52), labels = c("< -5", "-5 a -1", "0", "1 a 5", 
        "> 5"))

res.lm <- ggplot(setor.sf) + geom_sf(aes(fill = brks.res.lm), color = "black") + 
    ggtitle("Resíduos - LM") + scale_fill_discrete_sequential(palette = "Purple-Yellow", 
    c1 = 80, c2 = 30, l1 = 30, l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75", 
    drop = FALSE, name = "Taxa") + theme_void()

setor.sf$brks.res.car <- cut(dourados.car$residuals, include.lowest = TRUE, right = TRUE, 
    breaks = c(-14, -5, -1, 1, 5, 52), labels = c("< -5", "-5 a -1", "0", "1 a 5", 
        "> 5"))

res.car <- ggplot(setor.sf) + geom_sf(aes(fill = brks.res.car), color = "black") + 
    ggtitle("Resíduos - CAR") + scale_fill_discrete_sequential(palette = "Purple-Yellow", 
    c1 = 80, c2 = 30, l1 = 30, l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75", 
    drop = FALSE, name = "Taxa") + theme_void()


setor.sf$brks.res.gwr <- cut(results$residuos, include.lowest = TRUE, right = TRUE, 
    breaks = c(-14, -5, -1, 1, 5, 52), labels = c("< -5", "-5 a -1", "0", "1 a 5", 
        "> 5"))

res.gwr <- ggplot(setor.sf) + geom_sf(aes(fill = brks.res.gwr), color = "black") + 
    ggtitle("Resíduos - GWR") + scale_fill_discrete_sequential(palette = "Purple-Yellow", 
    c1 = 80, c2 = 30, l1 = 30, l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75", 
    drop = FALSE, name = "Taxa") + theme_void()


library(gridExtra)
grid.arrange(res.lm, res.car, res.gwr, ncol = 2)

Bibliografia sugerida

  • Druck, S.; Carvalho, M.S.; Câmara, G.; Monteiro, A.V.M. (eds). Análise Espacial de Dados Geográficos. Brasília, EMBRAPA, 2004.

  • Interactive Spatial Data Analysis by Trevor C. Bailey , Anthony C. Gatrell Routledge, 1995

  • Applied Spatial Statistics for Public Health Data; Lance A. Waller, Carol A. Gotway Wiley-Interscience 1St ed. 2004

  • Applied Spatial Data Analysis with R; Roger S. Bivand, Edzer Pebesma , Virgilio Gomez-Rubio Springer; Edição: 2nd ed. 2013

  • Geocomputation with R by Robin Lovelace, Jakub Nowosad and Jannes Muenchow. Geocomputation with R, 2021.

  • Spatial Statistics Workbook of Department of Criminology at the University of Manchester by Reka Solymosi and Juanjo Medina Crime Mapping in R

Obrigado !!!!!